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The dynamics of a trapped Bose-condensed gas at finite 
temperatures is described by a generalized Gross-Pitaevskii 
equation for the condensate order parameter and a semi- 
classical kinetic equation for the thermal cloud, solved us- 
ing TV-body simulations. The two components are coupled by 
mean fields as well as collisional processes that transfer atoms 
between the two. We use this scheme to investigate scissors 
modes in anisotropic traps as a function of temperature. Fre- 
quency shifts and damping rates of the condensate mode are 
extracted, and are found to be in good agreement with recent 
experiments. 
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The experimental observation of the scissors mode in a 
gaseous Bose-Einstein condensate (BEC) jjj provided a 
characteristic signature of superfluidity in this system . 
The scissors mode is excited by an angular displacement 
of the gas relative to an anisotropic confining potential. 
Above the BEC transition temperature, T c , the result- 
ing oscillation generally consists of the superposition of 
two modes, a high-frequency mode which reduces to an 
irrotational quadrupole mode in the limit of an isotropic 
trap, and a low-frequency mode which in the same limit 
corresponds to a pure rotation at zero frequency In 
contrast, a pure condensate exhibits only one oscillation 
frequency, indicating the irrotational, and therefore su- 
perfluid, nature of the gas. An important issue in super- 
fluid systems is the transition between these regimes with 
increasing temperature. This question has recently been 
addressed experimentally for trapped gases ||| , stimulat- 
ing the need for a consistent theoretical description of the 
observed behavior. 

It is well established that many properties of the con- 
densate at very low temperatures can be described by the 
Gross-Pitaevskii (GP) equation [||, which is a nonlinear 
Schrodinger equation for the condensate wavefunction, 
$(r,i). The equation treats the condensate as a classi- 
cal field and neglects quantum and thermal fluctuations. 
Consequently, the theory breaks down at higher temper- 
atures (T > 0.5T C ) where the non-condensed component 
of the cloud is significant. Including the thermal compo- 
nent in a consistent manner is a considerable challenge. 
Most calculations, such as those based on the Hartree- 
Fock-Bogoliubov (HFB) equations (|||, fail to capture 
the full collective dynamics of the thermal cloud, partic- 
ularly its back-action on the condensate. 

An approach which allows one to treat the dynamics 
of both components simultaneously was developed previ- 



ously [Q. The resulting equations of motion reduce to a 
generalized GP equation for the condensate 
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and a semiclassical Boltzmann kinetic equation for the 
thermal cloud 
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n c {r,t) — | 'I? (r, t) | 2 and n(r,t) are the condensate and 
non-condensate densities, respectively, and the mean field 
acting on the thermal atoms is given by U e s — V+2g\n c + 
h]. Apart from mean field effects, Eqs. (|l|) and (f|) are 
linked by the non-Hcrmitian source term, i?(r, t), which 
accounts for the transfer of atoms between the conden- 
sate and thermal cloud and is defined in terms of the 
Ci2[/] collision integral by 
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In arriving at this set of equations, several approxima- 
tions have been made. The use of a contact interatomic 
interaction, gS(r — r') (g = 4irh 2 a/m, where a is the 
atomic s-wave scattering length), is standard. More 
importantly, the thermal excitations are treated in the 
Hartree-Fock approximation. Furthermore, we invoke 
the Popov approximation, where the 'anomalous' density 
m is neglected. This is in part motivated by physical con- 
cerns: inclusion of m leads to unphysical low-momentum 
gaps in the energy spectrum, as well as infra-red and 
ultra-violet divergences. These inconsistencies can be re- 
moved either perturbatively or through the use of 
refined kinetic equations Q. The latter may provide a 
possible future extension of our work. 

In this paper we apply the above set of equations to 
the calculation of the scissors mode in an anisotropic har- 
monic trap [V(r) = m(uj 2 p 2 + uo 2 z 2 )/2] as a function 
of temperature for experimentally relevant parameters. 
The GP equation is solved using a FFT split-operator 
method ]lQ| , whereas the kinetic equation is solved by 
performing a classical simulation [ pL l|Jl 2f| in which the 
thermal phase-space density /(p,r,i) is represented by 
an iV-body system of discrete particles. Collisions be- 
tween atoms are naturally treated in this representation 



1 



by Monte-Carlo computation of the C22 and C\% colli- 
sion integrals. The complete dynamical problem becomes 
computationally feasible within this scheme, and has the 
additional advantage of providing an intuitive physical 
picture for the thermal gas dynamics. The damping rates 
and frequency shifts that we calculate for the scissors 
mode are in good agreement with experiment, indicating 
that our approach is valid over a wide range of tempera- 
tures. The simulations also yield the quadrupole response 
function of the system, which clearly illustrates the tran- 
sition between superfiuid and rigid-body behavior with 
increasing temperature. 

Our simulations involve propagation over a sequence of 
time steps, At. The trajectories of the thermal atoms are 
obtained by solving Newton's equations of motion with 
forces defined by the mean field potential U c g. To en- 
sure energy conservation the velocities and positions of 
the atoms are updated at each time-step using a second- 
order simplectic integrator jl3| . To define the mean field 
of a thermal cloud consisting of discrete particles, we first 
allocate atoms to grid-points using a cloud-in-cell ap- 
proach p3, and then convolve this density distribution 

2/2 

with a Gaussian function G(r) ~ e~ r > v . This effec- 
tively is a smoothening procedure, which is equivalent to 
assuming that interactions involving thermal atoms have 
a finite range rj. For consistency, the n c term appearing 
in U e s is also convolved. 

The Monte-Carlo computation of the C22 collision 
term, describing collisions between two non-condensed 
atoms that both remain in the thermal cloud, was dis- 
cussed in an earlier paper jla|. However, the C\i colli- 
sion term, which transfers atoms between the two com- 
ponents, was neglected. We now include this term, but 
as the method is somewhat involved we will detail it else- 
where |6). Briefly, one can write (|) as R = R in - R out , 
where R ln represents collisions between two thermal par- 
ticles that lead to absorption of one by the conden- 
sate, while R° ut refers to the inverse process. In terms 
of our Monte-Carlo simulations, one approximates the 
integrals by a sum over particles around a grid-point: 
R(r jkl ,t) = (V2n c )£i AP*, where AP t = P>» - P° ut - 
The term P™ (P° ut ) represents the probability that a 
given thermal atom with velocity Vj will collide dur- 
ing the time step, leading to a net transfer of parti- 
cles in (out) of the condensate. Both probabilities are 
proportional to n c and the Bose collision cross-section, 
a = 87ra 2 . They also depend upon the local condensate 
velocity v c and the phase-space densities of the final ther- 
mal states, randomly selected to satisfy momentum and 
energy conservation. 

In summary, the simulations consist of the following se- 
quence. In a given time step, the thermal particle phase- 
space coordinates are first updated as discussed above. 
The probabilities of C12 collisions are next calculated for 
each atom. A random number R € [0, 1] is chosen, and 
if |APj| > R, atoms are added (APj < 0) or removed 
(APj > 0) appropriately. The quantity APi is also ac- 



cumulated to define the dissipative term in ([!]) . We then 
treat the C22 collisions jl5|| . Finally, the GP equation is 
propagated; the dissipative term, R, leads to a continu- 
ous change in normalization of the wavefunction which is 
consistent with the discrete addition or removal of parti- 
cles from the thermal cloud. 

Experimentally the condensate is produced in a 
disk-shaped anisotropic trap (u> z ~ y/8ui p ), which is adi- 
abatically tilted to make an angle 0q with respect to its 
original orientation. The scissors mode is then excited by 
suddenly switching the trap to an angle — 9q, so that the 
condensate and thermal cloud oscillate about this new 
equilibrium position. We simulate this scenario by first 
finding the equilibrium condensate and non-condensate 
density profiles for a particular temperature T using a 
self-consistent semiclassical procedure [jl5). A sample 
of test particles is then chosen to simulate the phase 
space distribution, where to minimize statistical fluctu- 
ations, ten times the physical number of thermal atoms 
is used. Initial particle positions and momenta are ran- 
domly selected by a rejection method from a Bose dis- 
tribution f(p,r,t) = [z^ 1 exp(/?p 2 /2m) — 1] _1 , where 
[3 = l/k B T and z(r) = exp(-/3{[/ eff (r) - fi}) is the 
position-dependent fugacity (fx is the condensate chem- 
ical potential). The particle coordinates, as well as the 
condensate density, are then rotated through an angle 
28o about the y axis relative to the trap potential V(r). 

Starting with these initial conditions, quadrupole mo- 
ments can be calculated separately for the conden- 
sate (Q c (t) = fdrxzn c ) and thermal cloud (Q(t) = 

Ylj—i XiZi). To make contact with experiment we define 
a rotation angle for each component, 8 a (t) = Q a (t)/Q°, 
where Q° = (x 2 — z 2 )° a is the equilibrium quadrupole mo- 
ment of the a-th component [ |17|] . For a pure condensate 
(T = 0) consisting of N = 2 x 10 4 atoms, the quadrupole 
moment is found to oscillate with almost constant am- 
plitude (where small fluctuations arise from weak exci- 
tation of other condensate modes) at a single frequency 
of u! sc — 2.9886cj p , which is about 1.5% larger than the 
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nite number effects. Above T c the Bose gas oscillation 
exhibits two frequencies with approximately equal am- 
plitudes. Our simulations yield frequencies that are very 
close to those found experimentally and predicted ana- 
lytically: u>± — \lu p ±uj z \. In addition, the thermal cloud 
oscillation is weakly damped by C22 collisions, over a 
timescale which is similar for both modes and is of the 
order of several collision times. 

Below T c , our simulations describe the dynamics of 
both components. For most temperatures we find that 
the condensate and thermal cloud modes are essentially 
excited independently, indicating that the two compo- 
nents are only weakly coupled. Nevertheless, the conden- 
sate oscillation experiences significant damping from in- 
teractions with the thermal component. To quantify the 
damping rate and frequency of the condensate mode, the 
data is fit to a single exponentially-decaying sinusoidal 
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function in order to make contact with the experimental 
analysis H) . Fig. |l| shows results for the condensate mode 
as a function of temperature, for a fixed total number of 
atoms, N — 5 x 10 4 . To separate the effects of each term 
in the Boltzmann equation (^J) on the dynamics, simu- 
lations are performed with no collisions, with only C22 
collisions, and with both C12 and 622- At low tempera- 
tures (T < 0.6T®) we see that the damping is predomi- 
nately due to collisionless Landau damping, where mean 
field interactions between the condensate and thermal 
atoms transfer energy from the collective mode to single 
particle excitations. When C22 collisions are included, 
an increase in damping is observed. These collisions can 
only affect the condensate indirectly through their equili- 
brating effect on the nonequilibrium distribution of ther- 
mal atoms. The C12 term leads to additional damping 
at low T by further promoting the equilibration of the 
condensate and thermal cloud. This source of damping 
is comparable in magnitude to that of the C22 collisions 
but is small compared to Landau damping, in agreement 
with a previous study of collective modes under the as- 
sumption of a static thermal cloud jl8| . The effect of C12 
collisions can also be inferred by monitoring the number 
of condensate atoms as a function of time. Since the ini- 
tial angular displacement of the cloud increases its total 
energy, the gas must equilibrate to a higher final temper- 
ature with a corresponding reduction in the condensate 
number. This is seen in the simulations. 

At higher temperatures (T > 0.6T)?), collisional ef- 
fects increase in importance, but the situation becomes 
more complicated as Tj? is approached. The more mas- 
sive thermal cloud begins to drive the condensate at its 
own scissors mode frequencies, and as a result, a single 
damped sinusoid is a poor fit to the theoretical data in 
this regime. Effectively, the condensate oscillations are 
sustained over longer times and the damping appears to 
saturate. 

There are also observable effects on the condensate 
mode frequency. In the collisionless limit, we see from 
Fig. 1(a) that the frequency at first increases with in- 
creasing temperature as a result of the decreasing number 
of condensate atoms (i.e., the frequency shifts away from 
the TF limit). This trend is reversed by collisions, as 
well as at higher temperatures where the thermal cloud 
becomes significant. Since the condensate tends to drag 
part of the thermal cloud along with it, its effective iner- 
tia is increased. One would therefore expect a lowering 
of its normal mode frequency, as observed. 

Use of evaporative cooling in the actual experiments 
H meant that the total number of atoms varied with 
temperature from N ~ 2 x 10 4 at low T to N ~ 10 5 
close to Tj?. Since a small but significant dependence on 
number is found in our full simulations, we have taken 
this variation into account in our comparison with exper- 
iment. Our results in Fig. ^ are in very good agreement 
with experiment for T < 0.8T C ° but deviate from experi- 
ment in both frequency and damping at the three high- 
est temperatures. However the experimental error bars 



are particularly large for these points, which may reflect 
difficulties in extracting values of the damping rate and 
frequency from fits to the data over a limited timescale. 
We note that three damped sinusoids are needed to fit 
the theoretical data in this regime, as the condensate is 
strongly coupled to the thermal cloud. In addition, we 
do not see a condensate at the highest temperature point 
(T ~ which may indicate systematic errors in the 
experimental temperature scale or possibly limitations of 
the semiclassical approximation. 

A quantity of considerable interest is the quadrupole 
response function x"(uS), which physically describes the 
energy absorption of each mode under a harmonic per- 
turbation. Zambelli and Stringari JljJ demonstrated that 
this can be related to the moment of intertia of the sys- 
tem, O, by the expression 
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where O r i g id is the moment of inertia of the correspond- 
ing rigid body. For a sudden rotation of the trap, 
x"{w) cx u> %l{Q(uj)}, where Q(u>) is the Fourier trans- 
form of the time-dependent quadrupole moment. We ex- 
tract this quantity by fitting three damped sinusoids to 
the calculated total quadrupole moment. Fig. || shows 
results as a function of temperature. At lower temper- 
atures, the condensate mode is dominant and x"(w) ex ~ 
hibits a single peak. For very low T the thermal cloud is 
in fact strongly coupled to the condensate, and oscillates 
at the condensate frequency with a small phase shift that 
accounts for condensate damping. However, as the tem- 
perature is raised, the strength of the condensate mode 
diminishes, and the spectral density is dominated by the 
two thermal cloud modes at u)±, as would be expected 
for O ~ ©rigid p7[ ]. In contrast, the small T behav- 
ior is consistent with a superfluid moment of intertia, 
s f < ©rigid- In principle, one could calculate the mo- 
ment of inertia over the entire temperature range using 
(||). However, the damped nature of the modes leads to 
Lorentzian spectral densities for which the required fre- 
quency moments are undefined. Nevertheless, one can see 
qualitatively that @ s f < < ©rigid in the intermediate 
region. 

To summarize, we have simulated the scissors modes 
in a finite-temperature Bose gas using a coupled 
FFT/Monte-Carlo scheme. Our approach takes into ac- 
count fully the dynamical mean fields acting between the 
condensate and thermal cloud as well as all collisional 
processes that are physically relevant. We find very good 
agreement with experiment over a wide range of temper- 
atures, with only minor discrepancies near to T® where 
both theory and experiment are more difficult to ana- 
lyze quantitatively. These same methods can be used to 
study other collective modes, as well as problems such as 
condensate growth. 
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FIG. 1. Frequency (a) and damping rate (b) of the 
condensate scissors mode as a function of temperature for 
N — 5 x 10 4 total particles. Note that the temperature is 
normalized by the transition temperature of the ideal gas, 
T®. Each plot shows the effect of inclusion of each term in 
the kinetic equation (Q): the free-streaming operator alone 
(*); the free-streaming and C22 terms (•); and all terms (o). 
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FIG. 2. Frequency (a) and damping rate (b) of the scis- 
sors modes for a variable total number of atoms, intended to 
simulate the experimental parameters. The condensate mode 
is indicated by open (theory) and solid (experiment) circles. 
The open squares in (b) show the calculated average damping 
rate of the two thermal cloud modes, while the solid squares 
are the corresponding experimental values. 
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FIG. 3. Quadmpole response function (uj) f° r the sys- 
tem, as a function of frequency and temperature. The param- 
eters match those in Fig. ti. 
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